Multiple Invaded Consolidating Materials 
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We study a multiple invasion model to simulate corrosion or intrusion processes. Estimated values 
for the fractal dimension of the invaded region reveal that the critical exponents vary as function 
of the generation number G, i.e., with the number of times the invasion process takes place. The 
averaged mass M of the invaded region decreases with a power-law as a function of G, M ~ C 9 , 
where the exponent j3 w 0.6. We also find that the fractal dimension of the invaded cluster changes 
from di — 1.887±0.002 to d s = 1.217±0.005. This result confirms that the multiple invasion process 
follows a continuous transition from one universality class (NTIP) to another (optimal path). In 
addition, we report extensive numerical simulations that indicate that the mass distribution of 
avalanches P(S, L) has a power-law behavior and we find that the exponent r governing the power- 
law P(S, L) ~ S~ T changes continuously as a function of the parameter G. We propose a scaling 
law for the mass distribution of avalanches for different number of generations G. 

PACS numbers: 61.43.Gt, 05.40.-a, 64.60.Ak, 45.70.Ht 



I. INTRODUCTION 

The veins of gems and ores are often the product of a 
multiple intrusion of a reacting fluid into a porous soil in 
which dissolution and subsequent re-crystallization pro- 
cesses are the determining factor. Some examples like 
porphyry copper deposits 1] or olivine have been stud- 
ied in the literature and it is known that the surviving 
network of ore deposits has a fractal structure 0, Q that 
can be exploited for mineral exploration 5]. A similar 
situation can be found in vulcanology when magma is 
repeatedly injected through the same pathway, each time 
melting up again the most recent formations to find its 
way out |3 . 

The evolution of the pore structure after several 
invasion-frost-thaw events has been investigated numer- 
ically 0, and results indicate that the fractal dimension 
of invasion clusters varies with the number of invasion 
cycles. In this work, after invasion takes place, the struc- 
ture of the porous pathway is randomly healed. In a 
similar approach 8] , an optimized version of the multiple 
invasion percolation model was studied. Some topologi- 
cal aspects as the acceptance profile and the coordination 
number were investigated and compared to those of or- 
dinary invasion percolation. 

In the cases mentioned Q, 0, 0, 0, Q above and also in 
other cases of repeated invasions of corroding, dissolving 
or melting fluids into a strongly heterogeneous substrate, 
slowly consolidating matrix fractal patterns are created 
that reflect the history of the material. It is the aim 
of this paper to develop a model of multiple invasion in 
order to simulate how these patterns form and how their 
fractal dimension changes. In fact, we propose a complete 
theoretical framework based on scaling laws |(J . 

The theory of avalanche dynamics has been studied in 
a variety of contexts, for example in growth models, inter- 
face depinning and invasion percolation [jj. The forma- 
tion of fractal structures, diffusion with anomalous Hurst 



exponents and Levy flights, can all be related to the same 
underlying avalanche dynamics 0- Normally, the pres- 
ence of avalanches in the invasion process supposes un- 
changed porous media. In this work we also investigate 
the mass-distribution of avalanches and determine how 
the exponent that characterizes this distribution changes 
for different cycles of the invasion process. This paper is 
organized as follows. In Sec. II, we present the model to 
simulate the multiple invasion in consolidating medium. 
In Sec. Ill we show the results for the invaded cluster 
mass. The results and analysis of the numerical simu- 
lation for avalanche distribution are shown in Sec. IV., 
while the conclusions are presented in Sec. VI. 



II. MODEL 

In order to simulate the injection process we use the 
standard non-trapping invasion percolation (NTIP) |10| . 
In this model the invaded solid is considered to be very 
heterogeneous and the invading fluid can potentially en- 
ter anywhere along the interface. Here the consolidating 
medium is represented conveniently as a square network. 
The sites of the lattice can be viewed as the smallest units 
of constant strength and the randomness of the strength 
of the medium is incorporated by assigning random num- 
bers to sites. For simplicity, we consider the case in which 
dissolutions control the fluid invasion. 

On our heterogeneous medium we start by applying the 
standard invasion process of NTIP. For completeness the 
algorithm is described as follows. Initially, let us assign a 
random number, pi drawn from a uniform distribution in 
the interval [0, 1], to each site i of the lattice. We choose 
one site in the center of the lattice and occupy it. This 
site represents the injection point of the fluid and is the 
seed of the invading cluster. We look among the neigh- 
boring sites of this cluster (the growth sites) and choose 
the one which carries the smallest random number. This 
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FIG. 1: Typical cluster for different generations on a 256 x 256 
lattice, for p* = 1. The injection point is localized in the 
center of the lattice. In figures we show the results for G = 1 
and 100. 



site is then invaded and added to the cluster. Then we 
increase the list of sites that are eligible to be invaded. 
At each step of the invasion process, the perimeter of the 
nearest neighbors of the sites that form the invading clus- 
ter is investigated and the site that has the smallest pi is 
chosen. This procedure is repeated until the edge of the 
lattice is reached. At this point the simulation stops and 
the mass M (i.e. , the number of sites belonging to the 
invaded cluster) of the cluster is computed. The number 
of sites of the invaded cluster is very often considered as 
a time parameter. 

Now we present the new feature introduced to the stan- 
dard invasion percolation. After we finish the above de- 
scribed simulation in agreement with customary NTIP, 
the simulation is performed again starting every time at 
the same injection point. New random numbers chosen 
from a uniform distribution in the interval [0, 1] are as- 
signed to all sites belonging to the previously invaded 
cluster before a new invasion process starts. To all other 
sites, i.e., namely, those that are outside the cluster, we 
assign a random number homogeneously distributed in 
the interval [p* , 1] where p* is a number close to unity. 
Compared to the support used in the first generation 
where all sites can be invaded, the second generation 
appears substantially reduced, because it mostly corre- 
sponds to the cluster invaded in the first generation. In 
this way we generate again an invasion cluster for which 
p* = 1 is a subset of the previous one and so necessar- 
ily smaller. This procedure is repeated G times, where 
G is the number of generations. Standard invasion per- 
colation coincides with the case G = 1. At each new 
generation G, the sites of the previous invasion are re- 
invaded. 

Let us first consider the case p* — 1. In this situa- 
tion, the invaded cluster is after each time a subset of 
the previous one so that after a finite number of itera- 
tions the cluster does not decrease any longer. The num- 
ber of generations needed to reach a cluster whose mass 
remains unchanged depends on the size of the original 
lattice, because the number of possible available sites is 
proportional to the system size. Therefore, the satura- 
tion number is different for each lattice size. In order to 
illustrate these changes in the structure after each process 
of invasion, we show in Fig. ^ typical clusters generated 
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FIG. 2: The probability distribution Pi nv (p) of invaded sites 
for different generations G — 1 (circles), 2(squares) and 16 
(triangles), L = 512 and p* = 1. 



for a lattice of size L — 256 for four different generations 
G. Another important quantity is the probability distri- 
bution of pi of the invaded sites. In Fig.|2we present the 
normalized distribution Pi nv (p) for different generations 
G obtained from 1000 realizations of size L = 512. After 
the completion of the first invasion process, the distribu- 
tion expectedly displays a transition at p w p c , where p c 
is the critical site percolation point, p c = 0.59275 for a 
square lattice ^lj ■ The same behavior has been observed 
by numerical simulation in Refs. 0,0- For G = 2 the 
distribution Pi nv (p) becomes flat and the profile does not 
change any more as function of G. This happens because 
when G > 1 sites with larger pi are also invaded. 



III. CLUSTER MASS 

In our simulations we used the NTIP algorithm for 
square lattices of sizes L = 64, 128, 256, 512 and 1024. 
For each value of G, we perform simulations for 10000 
realizations and compute the mass Mq of the invaded 
cluster. In Fig. [31 we show the ratio Mq/Mq~i as a func- 
tion of the generation number G. For each size L, G s is 
defined as the number of generations at which the mass 
of the invaded cluster reaches a constant value, i.e., for 
which MgJMg s -i = 1. The results of our simulations 
shown in Fig. 0] for four values of the generation number 
G, indicate that the mass M has a power-law dependence 
on the size L, M ~ L d ° , where do is the fractal dimen- 
sion of the invaded cluster. The case G = 1 corresponds 
to the standard invasion percolation model. The value 
obtained from our simulations, d\ — 1.887 ± 0.002, is in 
good agreement with the current estimate d\ — 1.8959 
for NTIP [H El El 13. The results shown in Fig. H 
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FIG. 3: The evolution of the rate Ma/Ma-i as function of the 
logarithm of the generation number G for p* = 1. Here Mg 
is the mass at generation G for different sizes L — 64, 128, 256 
and 512. 
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follows a power-law behavior 

M ~ G" 3 . 



(1) 



To better analyze the data, we normalize the mass by the 
constant M\, which is the average mass of the invaded 
cluster at G = 1. Similar to some problems that involve 
growth surfaces ^(|, this process has two characteristic 
regimes: (i) power law evolution and (ii) saturation when 
G — > oo. To describe this behavior we propose the scaling 
relation 61 



M(G,L) = (G- Y„ 



Mi 



(2) 



where Nq is an off-set value for the generation number 
and a and z are scaling exponents. We assume that the 
scaling function f(x) has the form f{x) ~ x@ in the limit 
x <C 1 and f(x) = const when x 3> 1. Furthermore, 
a direct relation among exponents a, (3 and z can be 
obtained. We find M/M\ ~ C 3 for L » 1 and, since 
M/M\ ~ L Q in the saturation regime (G >> 1), we 
obtain that a = d s — d\ . 

In the crossover region, when the fractal dimension 
goes from d\ (G = 1) to d s (G = G s ) we have 



(G - N ) ~ L z . 

From these relations, we obtain that 

a 



(3) 



(4) 



and from the fact that the fractal dimension has reached 
the saturation value d s w 1.22, it gives a = —0.68. The 
inset of Fig.Elshows the data collapse obtained by rescal- 
ing M/Mi and G according to the scaling form Eq. J2J. 
In this case the best fit to the data gives (3 ~ 0.6. Sub- 
stituting into Eq. 01 we find z = 1.13. 



IV. AVALANCHE DISTRIBUTION 



FIG. 4: Log-log plot of the mass M of the invaded cluster 
versus the system size for different generation numbers G = 1 
(circles), 100 (squares), 500 (diamonds) and 3000 (triangles), 
and p* = 1. The straight lines are best fits to the data and 
their slopes are the fractal dimensions of the invaded clusters. 



indicate that by increasing the generation number the 
fractal dimension decreases continuously until it reaches 
a saturation value of d s = 1.217 ± 0.005 at G s . This 
value agrees with the fractal dimension of the optimal 
path in the strong disorder limit d opt — 1.22±0.01 [l5| . 
As shown in Fig. [3] for large system sizes we find that 
the average mass of the invaded cluster asymptotically 



It has been known since a long time that avalanches 
occur in invasion percolation and that these avalanches 
obey scaling relations related to percolation theory 
An avalanche occurs when a site j is invaded at a value 
Pj and then a series of sites i connected to this origi- 
nal site are sequentially invaded with pi < pj. It is also 
known that the system reaches a self-organized critical 
state characterized by avalanches of all sizes distributed 
according to a power law. In the case of NTIP, the ex- 
ponent corresponding to the power law behavior for the 
distribution P(S) of avalanche sizes S is r = 1.527 [l7|. 

In our simulation we found that the exponent corre- 
sponding to the case G — 1 is r = 1.46 ± 0.03. The 
expected value [17] is outside of our error bars, which 
we attribute to the fact that we have not reached the 
asymptotic limit because our systems are too small. 
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FIG. 5: Log- log plot of the average mass of the invaded cluster 
M normalized by the mass Mi of the first invaded cluster, 
against the number of generations G for different sizes L = 
64, 128, 256, 512 and 1024, and p* = 1. The inset shows the 
collapse following the scaling relation of Eq. ©• 



FIG. 6: The mass distribution of avalanches for different gen- 
eration numbers G = 2,4, 8, 16, 32, 64 and 128 for L = 512 
and p* = 1. The slopes of the straight lines follow power-laws 
with exponent r. The solid lines indicate the two limit cases 
G = 1 (lower) and G = 128 (upper). 



We performed simulations for different generations G 
on lattices of sizes L = 64, 128, 256, 512 and 1024, and 
calculated the size distribution of avalanches. In Fig.[5]we 
show P(S) for size L = 512 and G = 2,4, 8, 16, 32, 64 and 
128. It is clear from this figure that P(S) displays power- 
law behavior with exponent dependent of the number of 
generations G. The solid lines indicate the slopes in the 
two limit cases G = 1 (lower) and G = 128 (upper). 

In Fig. \7\ w e show how the exponent of the power- 
law t changes as function of the number of generations 
G. For large values of G the exponent converges to r = 
1. This value is the same found for the distribution of 
avalanches P(S) in the one-dimensional case |1S|. This 
is consistent with Fig. ^ for G = 100 where the avalanche 
process is limited to a thin path that is essentially an 
one-dimensional topology. 

In Figs.|HIa) andlHJb) we show the log- log plot of the 
distribution of avalanche sizes. It is clear from these fig- 
ures that P(S) displays a scaling region for intermediate 
avalanche sizes. In addition the scaling region is followed 
by a sudden cut-off that decays faster than exponential 
due to a finite size effect. The range of the power-law re- 
gion is proportional to the lattice size. As a consequence 
the biggest avalanches occur in the largest lattice. The 
position of the cutoff depends on G for fixed L. We pro- 
pose a scaling form for the mass distribution P(L,S), 
which accounts for finite size effects and power-law be- 
havior fijl 
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FIG. 7: Log-linear plot of the avalanche exponent r as func- 
tion of the generation number G, for L = 512 and p* = 1. 



where the function f(x) has a Gaussian form 

f(x) = exp[-x 2 }. (6) 

In practice, the appropriate parameters of the scaling 
function Eq. JSJ have been determined here through a 
nonlinear fitting procedure of the function 



P(S,L) = AoS-Texpi-iS/Ax) 2 } 



(7) 



(5) 



to the avalanche data. We observe that both the pref- 
actor Aq and the crossover amplitude A\ depend on the 
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FIG. 8: Log-log plot of the probability distribution of 
avalanches P(S, L) for various sizes L = 64, 128, 256, 512 and 
1024, p* = l.(a) G = 2 and (b) G = 128. 



system size. 

The solid line in Fig. corresponds to the best fit using 
Eq. Q for G — 2 and many different sizes L with r = 
1.37. The inset of Fig.[5]shows the power-law dependence 
of crossover amplitude on the system size, A\ a L 7 . The 
straight lines are the least-square fits to the data, with 
the slopes corresponding to the exponent 7 in Eq. JSJ for 
different generation numbers. 

In Fig. ^3 we plot the exponent 7 versus G, and see 
that the exponent has a monotonic behavior as function 
of the generation number. 

In Fig. ^2 we show the rescaled function P(S/L 1 ) for 
G = 16. The data collapse obtained confirms the valid- 
ity of the scaling form of Eq. @ . This confirms that the 
system is self-organized critical and the rescaled distri- 
bution shows the asymptotic scaling behavior of Eq. J7J . 
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FIG. 9: Log-log plot of the distribution P(S) for p* = 1 
and generation number G = 2 for L = 1024 (circles), 512 
(squares), 256 (up triangles), 128 (diamonds), and 64 (down 
triangles). The solid lines correspond to the scaling function 
y = A S~ T exp[-(S / A{) 2 ] with the parameter r = 1.37. The 
inset shows the log-log plot of the crossover amplitude Ai 
versus the system size L for G — 2 (circles), 32 (squares), 
128 (diamonds), and 256 (triangles). The lines are the least- 
square fits to the data and the slope is 7. 



V. RESULTS FOR p* / 1 

In the first part of this work we considered p* = 1. 
Now we present simulations for different p* very close to 
unity. In Fig. 1121 we show how the mass of the invaded 
cluster varies as function of the generation number G for 
a typical realization of the multiply invasion process. In 
the case p* = 0.9 the value of the mass shows strong 
fluctuations. If the probability to occupy sites outside 
of the previously invaded cluster is raised, the previous 
invaded region of the porous media is more likely invaded. 
To understand the qualitative behavior of the invaded 
cluster as function of the generation G we show in Figs. 1131 
and El typical clusters for two values p* = 0.9 and p* = 
0.999999, for five different generations G = 1,5,10,25 
and 50. For p* — 0.9, the cluster is more compact and 
sometimes changes the point where it reaches the border. 
When p* — 0.999999, the cluster becomes smaller at each 
generation. 

In order to be more quantitative we calculate the frac- 
tal dimension df. We measure the mass of the invaded 
cluster for different generations G for two different prob- 
abilities p* = 0.9 and 0.999999. Numerical simulations 
were carried out for 1000 realizations on lattice sizes 
L = 64, 128,256 and 512. In Figs. [T5I and [TBI we present 
log-log plots of the averaged mass of the invaded cluster 
versus the lattice size L. The linear fit to the data yields 
the fractal dimension df of the invaded cluster. In the 
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FIG. 12: Linear- log plot of the invaded mass for a typical 
FIG. 10: Log-linear plot of the exponent 7 versus the gener- realization as function of the generation number G, for L = 
ation number G. p* = 1. 512. From top to bottom, p* = 0.9, 0.9999, 0.999999. 
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FIG. 13: Typical configurations of invaded clusters at differ- 
ent generations G and L = 256. The random number, pi 
drawn from a uniform distribution of probabilities in the in- 
terval [p* , 1] for p* = 0.9. 



VI. CONCLUSIONS 



-6 -5 -4 -3 -2 -1 We have presented a comprehensive model to study 

log 10 S L a multiple invasion process. We have shown that the 

mass Mq of the invaded cluster decreases with the gen- 
eration number G. In addition, the fractal dimension 



FIG. 11: Log-log plot of the rescaled distribution of 
avalanches sizes P(S/L' y ) for generation number G — 16 
and different lattice sizes L = 64, 128, 256, 512 and 1204 and 
p* = l. 
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G=50 



case p* = 0.9, the fractal dimension is df = 1.90 ± 0.01 
for all generations. For p* = 0.999999 the fractal dimen- 
sion decreases when G increases. This implies that the 
fractal dimension of the invaded cluster has a behavior 
similar to the previously studied case in which p* = 1. 



FIG. 14: Typical cluster configurations for invaded clusters at 
different generations G and L = 256. The random number, 
Pi drawn from a uniform distribution of probabilities in the 
interval [p* , 1] for p* = 0.999999. 
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FIG. 15: Log-log plot the averaged mass M as function of the 
system size L for p* — 0.9 and G — 5 (circles), 10 (squares) 
and 50 (diamonds). The solid line with slope 1.90 ± 0.01 is 
the least-square fit to all data sets. 
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of the invaded cluster changes from d\ = 1.887 ± 0.002 
to d s = 1.217 ± 0.005 corresponding to G = 1 and 
G = G s , respectively. This result confirms that the 
multiple invasion process follows a continuous transition 
from one universality class (NTIP) to another (optimal 
path). We confirmed by extensive simulations that the 
invaded mass follows a power law M ~ G@ with an expo- 
nent j3 « 0.6. In addition the probability distribution of 
avalanches P(S,L,G) has been studied for different sys- 
tem sizes as function of the parameter G. We found that 
the mass distribution of avalanches follows a power law 
where the exponent r changes as function of the genera- 
tion number G. Based on this fact, we suggest that the 
avalanche process belongs to a different universality class 
for each G since no crossover scaling seems possible. Our 
results also indicate that this change in universality class 
occurs in a continuous way. Concerning the re-invasion 
of crystallizing, solidifying or healing fluids we conclude 
that only in the case in which the non-invaded part is not 
substantially damaged and the healed parts typically do 
not get much stronger than they were before the invasion, 
the multiple invasion process converges well to a differ- 
ent universality class, namely that of the optimal path 
13]. In the opposite case corresponding to p* ^ 1, the 
classical invasion percolation holds for all generations. 
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